function res = pfy (z, q)

   ## usage:  res = pfy (z, q)
   ##
   ## purify truncated normal

   global PAR

   persistent ipfy = 0 ;

   if (++ipfy >= PAR.pfy(1))
      res = z ;
      return ;
   endif 

   res = z ;
   n = rows(z) ;
   C = cov(z) ;
   z0 = repmat(norminv(q), n, 1) ;

   if (sum((z <= z0)(:)) / numel(z) < PAR.pfy(2)), return ; endif 

   [z flag] = gen1Z(C, z0) ;

   if !flag
      z = gen2Z(C, z0) ;
      if PAR.dbg, printf("pfy: using gen2Z.\n") ; endif
   endif 

   I = all(res <= z0, 2) ;   
   res(I) = z(I) ;

endfunction


function [res flag] = gen1Z (C, z0)

   ## usage:  [res flag] = gen1Z (C, z0)
   ##
   ## generate rows(z0) C randoms

   global PAR

   [nr nc] = size(z0) ;
   N = ceil(PAR.mem / (2*nc)) ;
   if (nc < PAR.pfy(4))
      z = mnormal_rnd(N, C) ;
   else
      if PAR.dbg, printf("pfy: gen1Z: Using truncated normals.\n") ; endif
      z = mnormal_rnd(N, C, z0(1,:)) ;
   endif 
   z0 = [z0 ; repmat(z0(1,:), N-nr, 1)] ;

   I = all(z <= z0, 2) ;
   res = z(I,:) ;
   iter = 0 ;
   # fid = fopen([tempdir "pfy_" num2str(getpid) ".log"], "wt") ;
   while (++iter < PAR.pfy(3) && rows(res) < nr)
      if (nc < PAR.pfy(4))
	 z = mnormal_rnd(N, C) ;
      else
	 z = mnormal_rnd(N, C, z0(1,:)) ;
      endif 
      I = all(z <= z0, 2) ;
      res = [res ; z(I,:)] ;
      # fprintf(fid, "%d  %10.2f%%\n", iter, 100*rows(res)/nr) ;
      # fflush(fid) ;
   endwhile 
   # fclose(fid) ;

   if (flag = (rows(res) >= nr))
      res = res(1:nr,:) ;
   endif

endfunction


function res = gen2Z (C, z0)

   ## usage:  res = gen2Z (C, z0)
   ##
   ## generate n C randoms (weak)

   [nr nc] = size(z0) ;
   z = mnormal_rnd(nr, C) ;

   res = z ;
   I = (res <= z0) ;
   res(I) = z(I) ;
   while any(!I(:))
      z = mnormal_rnd(nr, C) ;
      II = !I & (z <= z0) ;
      I = I | II ;
      res(II) = z(II) ;
   endwhile 

endfunction
